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A first order symmetric liyperbolic tetrad formulation of tlie Einstein equations developed by 
Estabrook and Wahlquist and put into a form suitable for numerical relativity by Buchman and 
Bardeen (the WEBB formulation) is adapted to explicit spherical symmetry and tested for accuracy 
and stability in the evolution of spherically symmetric black holes (the Schwarzschild geometry). 
The lapse and shift, which specify the evolution of the coordinates relative to the tetrad congruence, 
are reset at frequent time intervals to keep the constant-time hypersurfaces nearly orthogonal to the 
tetrad congruence and the spatial coordinate satisfying a kind of minimal rate of strain condition. 
By arranging through initial conditions that the constant-time hypersurfaces are asymptotically 
hyperbolic, we simplify the boundary value problem and improve stability of the evolution. Results 
are obtained for both tetrad gauges ("Nester" and "Lorentz") of the WEBB formalism using finite 
difference numerical methods. We are able to obtain stable unconstrained evolution with the Nester 
gauge for certain initial conditions, but not with the Lorentz gauge. 



I. INTRODUCTION 

An orthonormal tetrad approach to numerical relativ- 
ity has several attractive features. The metric is trivial, 
and most of the variables are coordinate scalars, which 
eliminates derivatives of the shift vector from most of 
the equations. There are only twenty-four connection 
coefficient variables in general, the Ricci rotation coeffi- 
cients, as opposed to the forty connection coefficients in a 
metric-based formulation. While one does have to evolve 
the tetrad vectors in place of the metric, one does not 
have to deal with nonlinearities in the equations associ- 
ated with the inverse metric. A number of tetrad and 
triad formulations for general relativity have been pro - 

posed QiiiiiBliiiaiiiiiiiiiiiiiaiaB, 

but they have not been as widely used in numerical rel- 
ativity as standard 3-1-1 metric-based formulations. This 
paper describes the numerical implementation and tests 
of one such scheme, the WEBB formulation j^, in spher- 
ically symmetric vacuum black hole spacetimes. See Es- 
tabrook's papers [isL Hgf for an in-depth analysis of the 
mathematical structure of these equations (without spec- 
ified gauge conditions). 

The WEBB scheme takes as its primary variables the 
tetrad connection coefficients, and incorporates one of 
two alternative dynamic gauge conditions to evolve the 
acceleration and angular velocity of the tetrads, either 
the Nester gauge or the Lorentz gauge 0- The evo- 
lution equations constitute a first-order symmetrizable 
hyperbolic system in which all variables propagate either 
along the light cone or along the tetrad congruence, with 
the tetrad gauge information propagating along the light 
cone. Both the Nester and Lorentz gauges do not in 
general preserve hypersurface orthogonality of the tetrad 
congruence. The coordinate evolution is controlled by a 
tetrad lapse function and shift vector which we do not 



evolve dynamically, but rather reset periodically to keep 
the constant-time hypersurfaces nearly (but not exactly) 
orthogonal to the congruence worldlines, and to maintain 
a minimal deformation condition on the spatial coordi- 
nates which is similar to the minimal strain condition of 
Smarr and York 21] in the 3-1-1 context. 

Spherically symmetric spacetimes, while in a sense 
trivial in that they do not allow any gravitational waves 
to be present, provide the challenge of maintaining a 
stable numerical evolution in the presence of an event 
horizon. Additionally, they require dealing with both 
an excision inner boundary and an outer boundary. In 
our code, the numerical grid extends from just inside the 
event horizon to around R — 20 M , where R is the cir- 
cumferential radius of a two-sphere. Our initial slices 
are constructed so that the congruence worldlines point 
out of the grid at both boundaries. This forces all of 
the eigenmodes to propagate out of the grid at the inner 
boundary, and the eigenmodes which travel along the 
congruence worldlines, as well those which travel along 
the outgoing light cones, to propagate out of the grid at 
the outer boundary. Even so, there are two eigenmodes 
(a "constraint" mode and a gauge mode) traveling at 
the speed of light into the grid at the outer boundary. 
Boundary conditions for the "constraint" mode are deter- 
mined according to constraint-preser ving boundary con- 
ditions m m ii m ie', "2?, "^s", "2?, "so; n m ni m, 

which insure that the information entering the numerical 
grid is consistent with the constraint equations. There 
is no physical constraint on the incoming gauge mode; 
we choose to keep its amplitude fixed as set in the initial 
conditions. 

This paper is one of a relatively small number of nu- 
merical tests of tetrad / triad formulations in vacuum gen- 
eral relativity (see IsM l3a 1371 ISSl l39| ) . We find that 
it is possible to achieve reasonable long term stability 
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evolving the spherically symmetric Schwarzschild geom- 
etry using the WEBB equations with the Nester gauge, 
but only for rather special initial conditions. Many other 
equation formulations and tetrad gauge conditions are 
possible in the context of orthonormal frames and remain 
to be explored. Tests in the limited context of spherical 
symmetry and one-dimension are far from sufficient to es- 
tablish the viability of a particular formulation, but can 
serve to establish lack of viability. 



II. VARIABLES 

The application of the three-dimensional WEBB 
formalism presented in to one-dimensional 

Schwarzschild black holes requires the construction 
of an orthonormal tetrad field well-behaved everywhere 
outside and in the vicinity of the event horizon. The 
most obvious basis vectors for a spherically symmetric 
spacetime consist of a timelike vector field Cq and a 
spacelike vector field e^, both orthogonal to the two- 
spheres generated by the symmetry, and two spacelike 
vector fields tangent to the two-spheres. The timelike 
vector field defines a timelike congruence of worldlines 
orthogonal to the two-spheres. The problem is that 
orthonormal vector fields tangent to a two-sphere cannot 
continuously be defined everywhere on the two-sphere. 
If these tangent vectors are chosen to be aligned with 
angular polar coordinates, and e^, they are degen- 
erate at the poles, 9 — and 6 = tt. These "spherical" 
spacelike basis vectors are unsuitable for direct use 
as the tetrad vectors in the WEBB tetrad formalism. 
There would be singularities in some of the connection 
coefRcients (Ricci rotation coefRcients) at the poles. 
Instead, we must define a "Cartesian" triad of spacelike 
vectors which is rotated from the spherical triad. A 
simple way to do this is to invoke the same rotation as 
a function of the polar angles that takes spherical to 
Cartesian basis vectors in flat space, namely: 
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62 = sin sin ef sin <j) cos 9 es + cos 4> e : 



(1) 



(2) 



However, the spatial triad vectors are not in general tan- 
gent to the constant-time hypersurface. As described in 
[l5j , we decompose the spatial triad vectors into a time- 
like component parallel to the congruence and a spacelike 
component tangent to the hypersurface: 



(4) 



The three- vector Ba is not a unit vector if Aa is not zero. 
The vectors and are tangent to the hypersurfaces, 
so Ag — = 0. Thus, Aa only has one degree of free- 
dom, in the ef direction. Af is the radial 3-velocity of 
a tetrad observer with respect to an observer at rest in 
the constant-time hypersurface. We can now write the 
spherical triad vectors as 

Bf = Af eo+B/ dr, eg = Bf dg, = B d^, (5) 

where r is the radial coordinate in the hypersurface. 
Since Bg and B ^ are unit vectors, they can be found 
from the metric. We find it convenient to define new 
symbols such that 
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with R the circumferential radius of the two-sphere. 

The directional derivative along the timelike vector of 
the tetrad can be related to coordinate derivatives by 
defining a tetrad "lapse" a and "shift" vector f3'^. The 
shift has only a radial component, so 



Do 



a 



(7) 



Unless the tetrad congruence is orthogonal to the 
constant-time hypersurfaces, the tetrad lapse and shift 
are different from the standard 3-1-1 lapse and shift. 

The directional derivatives along the spatial triad di- 
rections are 



Da 



d 



(8) 



for the Cartesian triad, which is related to the spherical 
triad by Eqs. O, (01, and (jJJ. The spherical directional 
derivatives are 



63 



(3) 



Df = Af Dn 



-X 



dr 



(9) 



While the Cartesian tetrad must be used to define the 
WEBB variables, once these are defined it is convenient 
to rotate back to the spherical triad to take explicit ad- 
vantage of the spherical symmetry. 

As in any Cauchy formulation for numerical relativity, 
the evolution of the spacetime is described by a sequence 
of spacelike hypersurfaces. Since the state of the sys- 
tem is specified on such a constant-time hypersurface, 
spatial derivatives must be evaluated at constant time. 



D, 



R 



de, D^ 



1 



R sine 



(10) 



The twenty-four connection coefficients are derived 
from the commutators of the Cartesian tetrad directional 
derivatives. In the WEBB formulation, a space-time split 
is made which groups the connection coefRcients into two 
3x3 dyadic matrices Nab = \ £bcd ^cda and Kab = ^boa, 
plus two spatial vectors, the acceleration a?, = Thoo and 



3 



the angular velocity (relative to Fermi- Walker transport) 
^f) = ^ £bcd Tdco of the tetrad frames. 

As a consequence of the underlying spherical symme- 
try, only the anti-symmetric part of the Nat is non-zero. 
The vector n with Cartesian components = ^ Satc N^c 
points in the radial direction, n = rif e^, and the con- 
straint equations arising from commutators of the spatial 
tetrad vectors reduce to 



1- Df R 
R ■ 



(11) 



The spherical symmetry also means that the Kab are 
symmetric, which is the condition that the tetrad con- 
gruence has zero vorticity and is orthogonal to some set 
of spacelike hypersurfaces. The Kab are the Cartesian 
components of the extrinsic curvature tensor of these hy- 
persurfaces (not necessarily the same as constant-time 
hypersurfaces). There are only two independent degrees 
of freedom in the extrinsic curvature, since the only non- 
zero components of the extrinsic curvature with respect 
to the spherical triad basis are 



spatial derivatives and source terms which are functions 
of the variables. The natural result of the WEBB for- 
malism is equations relating directional derivatives. Be- 
cause in general there are time derivatives hidden inside 
the spatial directional derivatives, we call the equations 
expressed in terms of the tetrad directional derivatives 
quasievolution equations (if they contain an explicit Dq) 
or quasiconstraint equations if they contain only spatial 
directional derivatives. (See the WEBB paper 15j for a 
complete discussion.) Because they are simpler, we first 
discuss the quasievolution and quasiconstraint equations. 

Since the Einstein equations are covariant under ro- 
tations of the tetrad, the quasievolution equations for 
Kii, Kt, and Uf, and the quasiconstraint equations for 
Kt and n^, can be derived directly from the Einstein 
equations in the e^, e^, basis. The quasievolution 
equations are 



Do Kr - Df 



Kr, Kan = Ki2 = K 



T- 



(12) 



Do Kt — Df rif — S_Kt, 



(16) 



(17) 



Finally, the angular velocity LOi, of a spherically sym- 
metric congruence is identically zero, and the accelera- 
tion can only point in the radial direction perpendicular 
to the congruence worldlines, so the only non-zero com- 
ponent relative to the spherical triad basis is af. The 
calculation from the Cartesian basis gives 



ai 



cosipsmt' flf, 02 = sin(psm6' Of, = cosW af. 



(13) 

The WEBB equations are greatly simplified if Kr, Kt, 
Uf, Uf, Bfi, Bt, and Af as used as variables instead of the 
Cartesian components. This results in a reduction of the 
total number of variables from thirty-six to seven. The 
lapse function and shift vector are not evolved dynami- 
cally, but are reset periodically to optimize the evolution 
of the coordinates. 

The spacetime metric is obtained by first calculating 
gfiiy _ ^jj. ^ rpj^g gfi.u ]-|2atrix is then inverted to 

give g^i,. The resulting metric is: 

= [-^2 ^ pr2 - aI) + 2 a Af] dr" 

+ 2e^ [aAf+ (3'' (1 - Aj)] dr dt 
+ e^^ (1 - A?) dr^ + i?2 dO^ + R^ sin^ 9 d(j>^. (14) 

The spatial metric of the constant-time hypersurface is 

df = e^^ (1 - A?) dr^ + R^ d0^ + R^ sin^ 9 dcf)^ . (15) 



III. TETRAD QUASIEVOLUTION AND 
QUASICONSTRAINT EQUATIONS 

The true evolution equations must be equations for 
partial time derivatives of the variables in terms of partial 



Do Uf — Df Kt = S_nf, 



(18) 



with 



S_Kr ^4~nj-Kl + K'? 



2 Uf 



S_Kt — h Kr Kt — Krp ~ Uf. — Uf Uf, 



s_ 



R 

Kt - Kr 
R 



af Kt + Uf {Kr~2 Kt). 



The momentum and energy quasiconstraint equations 
are, respectively, 

Df Kt = (^^n-KTHl-Ruf)^ ^^g^ 
it 



Dip Tl/f' — 



3 n? 2nf Kt 
~^ R 2~ 



{2Kr + Kt). (20) 



The gauge quasievolution equations for af are not co- 
variant under rotation, and must be derived from the 
Nester and Lorentz gauge conditions in a Cartesian basis 
(Eqs. (44) and (47) of 15]). The results are converted 
to our spherical basis variables. Both the Nester and 
Lorentz gauges give quasievolution equations for af of 
the form 



Dq af — Df Kr = S_af. 



(21) 



The source, S_af depends on the gauge. For the Nester 
gauge. 



S_af 



2 {Kr - Kt) 
R 



(22) 
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and for the Lorentz gauge, 
2 {Kr ~ Kt) 



S_af 



R 



2 {Kr nf + Kraf). (23) 



The Nester gauge quasiconstraint equation is trivial be- 
cause in spherical symmetry, the curl of a radial vector 
is zero. 



V. TRUE CONSTRAINT EQUATIONS 

To eliminate the time derivatives hidden in the quasi- 
constraint equations, we use the true evolution equations. 
It is convenient to take linear combinations of the result 
to get decoupled equations for Br dr Kt and Br dr Uf- 
We call these the true momentum and energy constraint 
equations. 



IV. TRUE EVOLUTION EQUATIONS AND 
THEIR HYPERBOLIC STRUCTURE 

In order to evolve Kr^ Kt, nf, and af numerically, the 
evolution equations must be expressed as partial deriva- 
tives along the coordinate directions r and t. To obtain 
these true evolution equations, Eqs. © and ((Tn|l are 
substituted into the quasievolution equations in Sec. IIIII 
Linear combinations of the quasievolution equations are 
taken to isolate the time derivative of each variable. The 
results can be lumped together in the form 



q+C^ Br a, q = S. 



(24) 



In this equation, 




l-Ai 



(Af I Q 

I Af Q Q 

Q Af I 

V I Af. 



and 



l-Ai 



S_Kr + Af S_af 
S_af + Af S_Kr 
S_Kt + Af S_nf 
S_nf + Af S_Kt 



The eigensystem of the characteristic matrix, , con- 
sists of four eigenmodes: two with eigenvalue 1/(1 — Af) 
and amplitudes af + Kr and Uf + Kt, and two with 
eigenvalue —1/(1 -I- Af) and amplitudes af ~ Kr and 
Uf — Kt- Decomposing the Dq operator into partial 
derivatives gives coordinates speeds 



si(r, t) 



S2ir, t) 



e^^ a 
1 + Af 



1-Af 



/3^ 



(25) 



(26) 



The two eigenmodes involving Kr and af , are "lon- 
gitudinal" modes, since they are constructed from com- 
ponents of the extrinsic curvature tensor and accelera- 
tion vector projected along Bf. The two involving Kt 
and Uf, are "constraint" modes, since Kt and rif are 
the variables which appear in the principal parts of the 
constraint equations. 



Br dr Kt - {Kr - Kt) 



(27) 



Br dr nf = - - Kr Kt - ^ 
+ Af Kt {af -I- Uf). 



(28) 



These equations are used in calculating initial conditions, 
boundary conditions, and as an accuracy check for the 
numerical evolution. 



VI. EVOLUTION AND CONSTRAINT 
EQUATIONS FOR Bf, Bg, AND Af 

Evolution equations for the triad vector components 
are also required. Recall there are only two indepen- 
dent components in spherical symmetry, Br = e^^ and 
Bt = l/R. Using Bt rather than i? as a variable in 
the numerics is motivated by a desire to maintain a form 
of the equations similar to what is necessary in a three- 
dimensional calculation, but also leads to a significant 
improvement in accuracy of the results for the Lorentz 
gauge. The evolution equations for Bf and Bx are 



{dt - Cp)Bf = - a Kr Bf 

Do Br = - {Kr 



drP"- 



a 



{dt-Ci3)B^^- aKTBg 

^ Dq Bj" — — Kj" Bj" ■ 



)Br, (29) 



(30) 



There is a constraint equation for Bt, 

Br dr Bt = - {Bt -Uf- Af Kt) Bt, (31) 

which is equivalent to Eq. (|11|) for the coordinate deriva- 
tive of R. This constraint is used to obtain Bt (and R) in 
the initial conditions, and as an additional check on the 
accuracy of the numerical evolution. The evolution equa- 
tion for Af follows from commuting Dq with the spatial 
directional derivatives, 



Dq Af — af — Kr Af — Br dr (In a) . 



(32) 
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VII. THE INITIAL VALUE PROBLEM (IVP) 

The initial value problem consists of finding values for 
Af, Kt, Uf, -Bti Bh, Kji, and af with which to begin the 
numerical evolution. We take the congruence orthogonal 
to the initial hypersurface, so Af = Q initially. KT^ nf, 
and R are obtained from Eqs. (|77|) . and (|^ with 

Af = 0. The first integral of the constraint equations, 



(1 -Rnf) 



1 



2M 



(33) 



where the constant of integration M is the Schwarzschild 
mass, makes numerical integration of Eq. H28(l for Uf un- 
necessary. Using Eq. (|33|l rather than Eq. (|28|1 to obtain 
Uf in the initial conditions has a substantial impact on 
the constraint errors early in the evolution, but after a 
few dynamical times leads to only a modest improvement 
in accuracy. 

The free initial data are Kn, af, and on the ini- 
tial hypersurface. The choice of Bfi is relatively trivial, 
since it sets the initial relationship of coordinate radius to 
proper radius and plays no role in the subsequent dynam- 
ics of the congruence or the geometry. With grid spacing 
uniform in coordinate radius, the initial choice of Br can 
affect numerical accuracy simply because it affects the 
relative grid resolution in different parts of the domain. 
Two particular prescriptions for Kr and Of are consid- 
ered: one which, along with appropriate choices for the 
lapse and shift, leads analytically to a time independent 
solution of the evolution equations (Time Independent 
IVP), and one which sets a uniform value for the trace of 
the extrinsic curvature on the initial hypersurface (Con- 
stant Mean Curvature IVP). 



A. Time Independent IVP 

The Schwarzschild spacetime has a time Killing vec- 
tor field (timelike outside the horizon, R > 2M), which 
means that coordinate systems can be found in which 
the metric is time independent. However, in the WEBB 
tetrad formulation the variables will be time indepen- 
dent only if the tetrad congruence is stationary as well 
as the (coordinate) metric. Reproducing an analytically 
time-independent solution is the simplest and most ba- 
sic test of a numerical code. Note that the standard 
Schwarzschild slicing giving rise to a static metric is not 
satisfactory for our purposes, since the normal tetrad 
congruence would be singular, with infinite acceleration, 
on the horizon. 

In a time-independent evolution, a congruence initially 
orthogonal to a constant-time hypersurface will be or- 
thogonal at all times, so Af will be zero at all times. 
This is consistent with Eq. (|32l) if the lapse is given by 



Br dr (In a) 



(34) 



Also, the coordinate time derivative of the geometrically 
defined curvature radius R must be zero, a condition 



which, by Eqs. (|30() and (|31|l . requires that the shift 
satisfy 



/?'■ = - 



a Br R Kt 
\ — R Uf 



(35) 



With dt Br = 0, Br dr {[3'' / Br) = -a Kr. 

Setting the partial time derivatives of Kr, and af to 
zero in their true evolution equations gives two simulta- 
neous equations for the proper radial derivatives of Kr 
and af. The equation for Kr is 



Br dr Kr = -{l-R Uf) 
R Kt S_Kr 



(1 — i? Uf) S_af 



{l-R Uf)^ - {R KtY 



(36) 



Note that S_af depends on the choice of tetrad gauge, 
ie. Nester versus Lorentz. 

Eq. ||2n|) is singular when {l-RufY -{R KtY = 0, or 
R — 2M (see Eq. (|33f) '). A solution regular on the horizon 
is obtained by requiring that the numerator of Eq. I|36|) 
vanish at i? = 2M, which implies a relation between the 
values of Kr and R Kt on the horizon (which we call 
Krh and Uq, respectively). For the Nester gauge, this 
relation is 



2M Krh 



1+4 [/q^- 8 |t/or 
-4Uo-8Uo \Uo\ 



and for the Lorentz gauge, it is 



2 M Krh 



-1 + 8 1^0 
4 Uo 



(37) 



(38) 



Uo is a free parameter, which must be negative for a 
black hole horizon. L'Hopital's rule can be used to find 
the starting value for Br dr Kr on the horizon. 

While one can integrate the equation for Br dr af 
obtained along with Eq. it is simpler to use the 

algebraic expression 



Uf = 



M/ R^ +RKt Kr 
{l-R Uf) ■ 



(39) 



Eq. is obtained by requiring that dt{R Kt) = 0, 

eliminating radial derivatives of R Kt and Uf using the 
constraint equations, eliminating the shift using Eq. I|35|l. 
and finally, simplifying with the help of Eq. (|33f) . With 
Qf given by Eq. (|39|l . we can determine the lapse every- 
where using Eq. H34|l . 

In summary, for a particular choice of Uo, one can 
obtain Kr, af, a, and everywhere so that the evo- 
lution of all the variables is time independent. The 
constant-time hypersurface generally does become sin- 
gular at some point inside the horizon, so the initial hy- 
persurface must be terminated at an excision boundary 
inside the horizon before af and/or Kr become too large. 
At large R, the generic behavior is that Kt and Kr ap- 
proach the same constant value. This constant is posi- 
tive if Uo is greater (less negative) than a certain critical 
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value, which is exactly —0.25 for the Nester gauge and 
about —0.29 for the Lorentz gauge. A positive at 
large R is highly desirable in dealing with outer boundary 
conditions, because then the shift at that outer boundary 
is negative and modes propagating along the congruence 
propagate out of, rather than into, the grid. Also, we 
find that expansion of the congruence along the radial 
direction {Kji > 0) everywhere is generally helpful in re- 
ducing growth rates of any unstable constraint-violating 
modes. 



B. Constant Mean Curvature (CMC) Slice 

An attractive slicing condition from the point of view 
of the conformal approach to the initial value prob- 
lem [43, ^3 Ql| is to impose a uniform value for the 
trace of the extrinsic curvature of the initial hypersur- 
face. This "Constant Mean Curvature" slicing, with 
(trace K = Kn + 2 Kt = Kq), allows testing of cases 
where the evolution of the hypersurfaces and the congru- 
ence is time dependent. In order to make the evolution 
of the congruence as "quiet" as possible, we choose the 
initial acceleration to satisfy the stationarity condition of 
Eq. (jSnj- Then, with Kr — Kq — 2Kt, the momentum 
and energy constraints can be integrated as differential 
equations for Kt and Uf, starting at the horizon, as de- 
scribed in the previous section. Uq (the value of R Kt on 
the horizon) is a free parameter, along with the choice of 
Kq. The value of rif on the horizon is fixed by Eq. 

Together with a shift given by Eq. (|35|) and a lapse 
given by Eq. H34|l . this implementation of the CMC ini- 
tial condition means that the time derivatives of Kt and 
Uf vanish on the initial hypersurface, though they will not 
stay zero. A critical test of the Nester and Lorentz tetrad 
gauge conditions is whether the hypersurfaces and con- 
gruence will subsequently evolve toward or away from a 
time-independent solution. Of course, the answer to this 
question also depends on how the coordinates evolve, as 
described in the next section. 



VIII. RESETTING COORDINATE 
CONDITIONS 

Our system of equations is symmetrizable hyperbolic 
for any fixed choice of lapse and shift. Accordingly, dur- 
ing a given time step, the lapse and shift are held fixed at 
the values they have at the start of the time step. While 
it is essential to accuracy and stability that the system 
be hyperbolic during the time steps 25], it is not critical 
that the overall evolution be hyperbolic, only well-posed 
fisf. In other words, the lapse and the shift can be reset 
at fixed time intervals according to conditions which do 
not necessarily preserve the hyperbolicity of the system. 
This gives a wide range of choices for resetting conditions. 
Keeping a fixed lapse and shift for many dynamical times 
is likely to give rise to coordinate singularities. 



The lapse determines the evolution of the constant- 
time hypersurfaces. The hyperbolic system breaks down 
if Af. = 1, which signifies that the constant-time hyper- 
surface has become null. We reset the lapse to keep Af 
small, which is accomplished by making the reset lapse 
satisfy Eq. (|34|l . Then the time derivative of Af will be 
small by Eq. (|32|l as long as Af is small. Since Af = 
initially, it does in fact stay very small if the lapse is re- 
set at small time intervals. The constant of integration 
in solving Eq. H34I) is conveniently chosen to keep the 
lapse constant at the outer edge of the grid, but makes 
no practical difference, since it is just a uniform rescaling 
of the time coordinate. Choosing the time step accord- 
ing to a Courant condition automatically compensates 
for any such rescaling. 

The shift controls the evolution of the spatial coordi- 
nates. We work with a grid at fixed values of the ra- 
dial coordinate. The evolution of the radial coordinate 
should be managed to keep the grid from being sucked 
up by the black hole, while keeping the inner edge just 
inside the event horizon, without any excessive stretch- 
ing or compression of the grid relative to the physical 
curvature radius R. If made possible by the evolution of 
the congruence and the constant-time hypersurfaces, we 
want our variables to approach a stationary final state. 
One option would be to choose the shift to make the 
partial time derivative of the curvature radius R zero ev- 
erywhere at each resetting. From Eqs. (|30|l and H31|) this 
condition is 



13'' = - 



a BrRKt 



l-Ruf- Af RKt 



(40) 



However, it seems more desirable to use a condition which 
is not so tied to the special circumstances of spherical 
symmetry. 

The minimal strain condition, as introduced by Smarr 
and York |2lj for 3+1 formulations, minimizes an integral 
of the square of the time derivative of the spatial metric 
over the hypersurface. 



h 



dhi-j dhki 



dt dt 



(41) 



with respect to variations in the shift. In /i, hij is the 
spatial metric and h is the determinant of this metric. 
The metric time derivatives depend on the shift through 



dt 



-2 a K,. 



-'(3 ^ij ; 



(42) 



where Kij is the extrinsic curvature of the hypersurface. 

In the context of the WEBB equations, the vectors Ba 
carry the information corresponding to the spatial met- 
ric in a 3+1 formalism. These vectors are not unit vec- 
tors with respect to the true spatial metric if the tetrad 
congruence is not orthogonal to the constant-time hyper- 
surface. As projections of the spatial triad vectors into 
the constant-time hypersurface, they can be thought of 
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as representing unit spatial displacements in the local 
3-space orthogonal to the tetrad congruence, which are 
Lorentz-contracted when measured in the hypersurface 
frame. We find it convenient to pose a minimal deforma- 
tion condition on the components Bj^ , which in general 
is not the same as the Smarr-York condition. 

Inverting the Bj' matrix gives the matrix B\ of com- 
ponents of one-forms dual to the Ba vectors. A spatial 
metric with respect to which the Ba are unit vectors 
is hij = B"^ . We use this metric, rather than the 
true spatial metric, in our minimal deformation condi- 
tion, minimizing with respect to variations in the shift 
vector the action 



dB,, 



dt dt 



(43) 



where h is the determinant of hij . In the current applica- 
tion, with Af kept very small by resetting the lapse and 
a diagonal Bj', I2 is very nearly equivalent to Ii. 
The dependence on the shift is through 



and 



with 



dtBn drBB. „ „ „r 
—S— = P —5 a Kr - dr (3"^ = y 



= p ct Kt- 



Bx 



Bx 



dr 47r 



Bu 



dtBn 
Br 



The Euler-Lagrange equation is 



dr V 



,drBj 
Bt 



drBj 

Bt 



dtBj 
Bt 



a Ki 



(44) 



(45) 



(46) 



(47) 



which, along with Eq. (|44|1 giving drl3^ in terms of y, can 
be integrated to find (3^ and dr(3^ everywhere. Note that 
drBT can be evaluated from the constraint equation (|31|l . 
Strict minimization requires y = at each boundary, but 
in order to keep the inner edge of the grid just inside 
the horizon and the outer edge at a fixed R we instead 
impose Eq. (|in|l as the shift boundary condition at both 
edges of the grid. After each recalculation the values of 
P"^ and drP"^ at each grid point are stored and kept fixed 
until the next resetting of the lapse and shift. 



IX. BOUNDARY CONDITIONS 

We now discuss boundary conditions for the evolution 
equations. The inner boundary of the grid is an excision 
boundary, kept at a roughly constant R < 2M, inside the 
event horizon, by the shift condition of the last section. 
All characteristics of the evolution equations propagating 
at less than or equal to light speed (which is all of them 



in the WEBB framework) are then outgoing relative to 
the grid at the inner boundary, since inside the horizon 
the entire future light cone is toward decreasing R. All 
information needed to update the variables at the inner 
boundary can be obtained from the boundary values or 
upwind differencing. 

The outer boundary is a different story. We do not 
try to extend the outer boundary all the way to infi- 
nite R, though this in principle is possible. Since our 
asymptotically hyperbolic hypersurfaces are asymptoti- 
cally null, the variables should be regular functions of 
Bt = i/R, and the radial coordinate could be scaled 
to be linear in Bt approaching null infinity. We termi- 
nate our initial grid at a value of R the order of 20 M, 
where variables such as Uf, Kt, and Kr have reached 
nearly constant values, and keep the outer boundary at 
approximately the same R with our shift condition. This 
allows the grid to be roughly evenly spaced in proper 
radius, with adequate resolution near the horizon and a 
reasonable total number of grid points. As long as Kt 
is positive near the outer boundary, the shift is nega- 
tive there (see Eq. H40|l 'l. and the modes propagating 
along the hypersurface normal are outward relative to 
the boundary. This is critically important, since there is 
no good way of specifying incoming boundary conditions 
for these modes. Still, the two modes with speeds S2 are 
incoming at the boundary and do require boundary con- 
ditions. The amplitudes of these modes are {af + Kr) 
and {uf + Kt). Hence, we refer to them as the incom- 
ing "longitudinal" and "constraint" modes, respectively. 
The "constraint" modes consist of the variables whose 
radial derivatives appear in the principal parts of the en- 
ergy and momentum constraint equations, and have the 
potential of propagating constraint errors in from the 
boundary (see '25| for a fuller explanation). The con- 
straints cannot be satisfied by just setting the incoming 
constraint mode amplitude to zero. A number of differ- 
ent constraint-preserving boundary conditions have been 
proposed 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, slIH. 
What we do is, first, iteratively correct the values of the 
constraint variables at the boundary physical points so 
that the energy and momentum constraints evaluated at 
half a step inside the grid are satisfied. Then all vari- 
ables are extrapolated to the ghost point just outside the 
boundary, and the energy and momentum constraints are 
again applied iteratively halfway between the boundary 
physical point and the ghost point to adjust the extrap- 
olated values of the constraint variables. 

For stationary initial conditions, the longitudinal vari- 
ables, Kr and a,-, are kept equal to their initial values 
at the ghost point, to give boundary conditions for the 
incoming longitudinal mode, Kr + af. For CMC ini- 
tial conditions, the ghost point longitudinal variables can 
either be extrapolated quadratically or cubically to the 
ghost point, which could be dangerous for stability, or 
again maintained at their initial values. The incoming 
longitudinal mode amplitude is a purely gauge quantity, 
so the only real concern is not to do something which 
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leads to a gauge instability at the boundary. 



X. NUMERICAL METHODS 

The numerical methods used in the codes are designed 
to be second order accurate in both space and time. 
While we have experimented with various methods, the 
most accurate results, and the ones presented in this 
paper, were obtained with a second order Strang split 
method. During each full time step, the source terms 
alone are integrated for a half time step at each grid point 
using second-order Runge-Kutta, followed by a full wave 
propagation time step with the updated variables solv- 
ing the equations with source terms omitted, and then 
another half time step of source term integration. In the 
wave propagation time step the equations are reshuffled 
into decoupled advection equations for the eigenmode 
amplitudes, of the form 



dty + v dry = 0, 



(48) 



where v is the wave speed. Using characteristic tracing, 
the solution of the advection equation for the updated y 
at the J*'' grid point is just the initial value of y at the po- 
sition ri — v At, where v is the average advection velocity 
along the characteristic. This average speed is approxi- 
mated by advancing v to halfway through the time step 
by a first-order accurate Euler method at all grid points 
before doing the wave propagation (only necessary for 
the modes whose wave speed depends on Bn), and then 
interpolating 1/v to the midpoint in r along the charac- 
teristic during the wave propagation, so at the i*'* grid 
point 



(49) 



with upper and lower signs for {vi > 0) and {vi < 0), 
respectively. Standard second-order methods interpo- 
late y to the characteristic location using an upwind 
first-order difference and a second-order difference ei- 
ther centered at i (Lax-Wendroff) or centered at the up- 
wind grid point (Beam- Warming) . Except at the bound- 
aries, we use a linear combination of the two second- 
order differences which does the interpolation to third- 
order accuracy, weighting the centered second-order dif- 
ference by (2 — |w|)/3 and the upwind second-order dif- 
ference by (1 -I- I^D/S. While the overall method is not 
third-order accurate, this does significantly reduce the 
errors compared with pure Lax-Wendroff. At the bound- 
aries either pure Lax-Wendroff or pure Beam- Warming is 
used as necessary to minimize use of ghost cells, Beam- 
Warming for outgoing modes and Lax-Wendroff for in- 
coming modes. Note that the upwind point changes as 
V changes sign, which means a small but sudden change 
in the interpolated upwind second-order difference and 
its contribution to the updated y. The effect of this is 
noticeable in the constraint error plots we present in Sec. 
IXII but has no impact on stability. 



The integration of the spatial ordinary differential 
equations during the initial condition and boundary con- 
dition subroutines is implemented with a simple second- 
order Predictor-Corrector scheme iterated until the 
differences between the predicted and corrected values 
are all near machine precision. 



XI. RESULTS 

All our numerical results are obtained on a uniform 
grid ranging from —0.16 < r < 9.84, where r is the co- 
ordinate radius, with r = at the initial location of the 
event horizon. The initial Bji = 1, so that the coordi- 
nate radius starts out equaling the proper radius. The 
natural logarithm of the lapse (In a), the shift, {(3^), and 
their derivatives are reset as described in Sec. IVIIII For 
convergence studies, these quantities are reset at a con- 
stant coordinate time interval as the grid resolution is 
increased and the time step is decreased in accordance 
with the Courant condition. Note that all values plotted 
are in units with 2 M =1. 



A. Time 

When presenting results, it is desirable to relate the 
coordinate time (t) of the numerical evolution to a phys- 
ically meaningful time. The natural choice for a physical 
time is the proper time of an observer at infinite distance 
from the black hole, i.e., the Schwarzschild time coor- 
dinate ts- Our spacelike hypersurfaces are quite differ- 
ent from the constant-time surfaces in the Schwarzschild 
metric, but we can still relate changes in our time coor- 
dinate to changes in Schwarzschild time. 

The derivation is outlined assuming Af = 0, since Af 
is kept close to zero with our lapse resetting condition. 
The proper time for a displacement perpendicular to the 
two-spheres as obtained from the Schwarzschild metric is 



dr 



1 - 



2 M 



R 



dtl- 



2 M 
R 



(50) 



from which, with dt the change in our coordinate time, 



dt 



1 - 



2 M 

R. 



1 - 



2 M 
R 



dRV 



(51) 

For a displacement at constant coordinate radius r the 
metric of Eq.((Ti|) gives {dr/dtf = (a^ - {P^- e^)^). At 
both edges of the grid our boundary condition on the 
minimal deformation shift (Eq. I40|) gives dt R = 0, and 
simplifying with the help of Eq. I|33|l we arrive at 



dts 



\ — R Tlf 



dt 



(52) 
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for the relation between between Schwarzschild time and 
coordinate time there. We can normaUze the lapse at 
each resetting so that a — {1 — R Uf) at the outer- 
most grid point, so that the change in Schwarzschild time 
equals the change in coordinate time at the outer edge 
of the grid to reasonable accuracy. For stationary initial 
conditions, this also makes dts = dt hold at all interior 
grid points. 



B. Stationary IVP 

Testing the code using a stationary solution as the ini- 
tial condition is the first step toward determining the 
viability of the code. Ideally, as the numerical evolution 
proceeds, all the variables stay constant. This does not 
occur in practice, however, because numerical errors act 
as perturbations to the stationary solution. Even with 
an accurate and stable numerical method, these pertur- 
bations may excite analytic unstable modes of the evolu- 
tion equations, some of which may be constraint- violating 
and some of which may be purely gauge instabilities. 

Figs, n and |21 show solutions of the stationary initial 
value problem for the Nester and Lorentz gauges, respec- 
tively. The solutions are one-parameter families charac- 
terized by L/qj the value of R Kt on the horizon. We 
only consider cases where the radial component of the 
extrinsic curvature Kfi is positive (expansion of the nor- 
mals radially) at large r, since only then do modes prop- 
agating along the hypersurface normal propagate off of, 
rather than onto, the grid at the outer boundary. 

For the Nester gauge, only the inner portion of the 
spatial domain is shown, from just inside the horizon at 
r = —0.16 to r = 4.0, since the asymptotic behavior of 
most of the variables is already apparent by r = 4.0. 
The outer edge of the grid is always at r = 9.84. At the 
critical value Uq = —0.25 all of the variables plotted ap- 
proach zero as i? — > cxD, and the constant-time hypersur- 
face asymptotically approaches a Schwarzschild constant- 
time hypersurface. For Uq less negative than the critical 
value, both and Kt approach the same limiting (pos- 
itive) value, with trace K 0.15 for Uq = —0.18 and 
trace K 0.24 for Uo = -0.14. The gradients in Kr 
and Qf become steeper in the vicinity of the event hori- 
zon at r = as C/q becomes less negative. Note that the 
first-integral constraint of Eq. (|33|l requires a change of 
sign of rif from positive to negative going outward when 
Kt approaches a non-zero value at large R. This occurs 
outside the range of these plots, but not too far outside 
for Uo = -0.14. 

For the Lorentz gauge, we show the entire range of the 
grid, from just inside the horizon at r = —0.16 to the 
outer edge at r = 9.84. The critical value of Uq at which 
the asymptotic values of Kt and Kr change sign is about 
—0.29. The asymptotic behavior at large R when Uq is 
less negative than the critical value is not well-behaved. 
Both Kfj and Kt get steadily larger, eventually in a run- 
away indicating that the stationary condition is forcing 



the hypersurface to become null. The beginning of this 
runaway is apparent in the plot for Uq — —0.25, but is 
well outside the edge of the grid for Uq ~ —0.27. This 
indicates a pathology in the Lorentz gauge which does 
not allow an approach to a stationary solution over the 
entire spacetime, but may not cause problems when solu- 
tions are evolved over a limited domain. The pathology 
is caused by the extra terms in the Lorentz gauge source 
for the af evolution equation (Eq. . Once rif becomes 
negative the —2 Kr Uf term becomes positive and exerts 
positive feedback on Kr through Eq. 

To test the accuracy and stability of the time evolution 
with stationary initial conditions we can not only check 
how accurately the constraint equations of Sections Ivl 
and IVII are satisfied, but also take advantage of the fact 
that any change in any of the variables with time can 
only be due to some error in the calculation. The growth 
of some the constraint errors with time for both Nester 
and Lorentz gauges is shown in Fig. |3| The mass con- 
straint error is the fractional difference between the mass 
as calculated at each grid point and each time from Eq. 

and the mass used in constructing the initial con- 
ditions (M = 1/2 for our choice of units). The energy 
constraint error is the difference of the left and right-hand 
sides of Eq. H28|) , with the derivative of Uf evaluated nu- 
merically. The momentum constraint error (not plotted) 
from Eq. (|27|l is similar in magnitude to the energy con- 
straint error. The energy constraint errors are weighted 
by a factor of R, which serves to make the errors more 
representative of fractional errors in individual terms of 
the energy constraint and also to make the amplitudes of 
the curves more uniform with radius. We multiply the 
actual errors by a factor of four because when testing con- 
vergence, we compare these curves with those from runs 
with double the point spacing. Assuming the quadratic 
convergence expected from our second-order accurate nu- 
merical method, the errors quadruple when the spacing 
between grid points is doubled. 

The first two of the three plots show the results for the 
Nester gauge with Uo = —0.14, up to a time of 120 M. 
The errors do grow somewhat with time, but the rate of 
growth, after an initial oscillation in the mass constraint 
error, is clearly decreasing. The errors stop growing al- 
most completely at later times, as shown in Fig. 0]for the 
energy constraint errors up to i = 400 M. The Lorentz 
gauge results for the energy constraint errors, the third 
plot in Fig. 121 clearly show exponential growth with time. 
The growth rate is not sensitive to numerical errors, and 
seems to be due to an analytically unstable constraint- 
violating mode of the evolution equations. 

The apparent stability of the time evolution for the 
Nester gauge is encouraging, but it is also important to 
demonstrate convergence of the errors as the grid spacing 
is reduced. Fig. [Sldoes this for the Nester gauge evolution 
at t = 40 M. We compare results for dr = 0.02, dr = 
0.01, and dr — 0.005. The errors are multiplied by the 
appropriate factors to make the curves lie on top of each 
other assuming perfect second-order convergence. While 
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there are some deviations close to the inner edge of the 
grid, the convergence is quite good everywhere else. At 
later times the convergence does gradually deteriorate, 
but that is not surprising after a very large number of 
dynamical times. 

The good results for the Nester gauge with Uq = —0.14 
unfortunately deteriorate rather rapidly as Uq is varied 
away from this value. For Uq ~ —0.19, for instance, 
the errors grow much more rapidly initially, and then 
start oscillating with slowly increasing amplitude for t > 
100 M or so. Convergence starts breaking down around 
t — 60M. The errors at late times are at least a couple 
of order of magnitudes greater than they are for Uq — 
—0.14. There is no real indication of exponential growth, 
but accuracy is not very good even for 1000 points in the 
grid. 

The apparent need to fine-tune the initial conditions to 
get errors under good control with the Nester gauge sug- 
gests that the present equation formulation does not give 
the kind of robust accuracy and stability that is necessary 
for dealing with physically more interesting problems of 
black hole dynamics and mergers. The Lorentz gauge, 
with its exponentially growing instabilities, seems clearly 
unsuitable for dealing with black hole problems. 



C. Constant Mean Curvature (CMC) IVP 

A more demanding test of the code is to start from 
initial conditions which do not correspond to a station- 
ary solution, and see if the gauge conditions allow relax- 
ation toward a stationary solution. Certainly one only 
expects to be able to evolve black hole spacetimes ac- 
curately for many dynamical times if this is the case. 
The obvious first choice, consistent with the conformal 
app roach to the initial value problem pioneered by York 
|40l I41I , is to require that the trace of the extrinsic 
curvature be uniform on the initial hypersurface, which 
is what we mean by constant mean curvature. This is 
the condition for separability of the momentum and en- 
ergy constraint equations in the conformal approach, and 
has been assumed in much of the work on initial value 
problem over the last few decades. The most common as- 
sumption has been that the initial hypersurface is max- 
imal, trace K ~ 0, but hyperbolic hypersurfaces, with 
trace K > 0, are more desirable from our point of view, 
since they simplify boundary conditions at the outer edge 
of the grid and seem to be necessary for the stability of 
our evolution equations. 

The value of the constant mean curvature, A'o, is a free 
parameter, as is the value of Uq (i? Kt on the horizon). 
While the initial acceleration of the tetrad congruence 
is in principle a free function, we fix it according to the 
stationary condition of Eq. (|39|l . This gives the time 
evolution a quiet start, in that the initial partial time 
derivative of R Kt is zero, but nothing more than that. 

Fig. shows the Nester gauge time evolution of af 
for three different choices of Kq and Uq. With Kq — 0.2 



and Uq = —0.14, the solution approaches a stationary 
solution, with the value of R Kt on the event horizon 
equal to -0.15 at t ^ 200 M. With Kq = 0.2 and Uq = 
—0.18, the solution also approaches a stationary solution, 
with the value of R Kt on the horizon equal to —0.16 at 
t = 200 M. However, with Kq = 0.4 and Uq = -0.14, 
the solution becomes pathological due to an increasingly 
sharp gradient of the acceleration somewhat outside the 
horizon, and the code crashes. 

Similar initial conditions with the Lorentz gauge do 
not lead to any pathologies in the acceleration, but in 
all cases roughly exponential growth of constraint errors 
prevents meaningful continuation of the solutions much 
beyond a time of 40 M. 



XII. DISCUSSION 

The WEBB tetrad formulation for numerical relativ- 
ity has been implemented and tested in spherically sym- 
metric black hole spacetimes, with two types of tetrad 
gauge conditions, which we call the Nester and Lorentz 
gauge conditions. While there is freedom in the choice of 
the initial velocity and acceleration of the tetrad congru- 
ence, subsequent evolution of the tetrad frames follows 
uniquely from the tetrad gauge conditions. The coordi- 
nate evolution is determined by elliptic equations for a 
tetrad lapse and shift in such a way that hypersurfaces 
of constant time stay perpendicular to the tetrad con- 
gruence, and the spatial coordinates evolve according to 
a minimal deformation condition, with boundary condi- 
tions which keep both the inner and outer edges of the 
grid at roughly constant Schwarzschild curvature radius 
R. While we pay particular attention to initial conditions 
consistent with a stationary evolution, we also consider 
non-stationary evolution from constant mean curvature 
initial hypersurfaces. Some representative results have 
been presented, but we have experimented with a variety 
of other parameter choices, numerical methods, etc. We 
have done sufficient testing to convince ourselves that the 
instabilities we find are not due to problems with numer- 
ical methods, but represent genuine analytic instabilities 
of the evolution equations, with numerical errors only 
playing a role of seeding the instabilities. 

There are several potential factors which can influence 
the existence and severity of instabilities, such as (1) the 
choice of tetrad gauge, (2) the equation formulation, in- 
cluding choices of variables and the possible addition of 
constraint equations to the evolution equations in vari- 
ous combinations, (3) initial conditions and how the co- 
ordinates are evolved, and (4) boundary conditions. Of 
course, all of these factors interact in various ways and 
it is not always clear just what is behind a particular 
instability. 

The equation formulation, gauge conditions, and the 
choice of variables determine the hyperbolic structure 
of the evolution system. The WEBB framework makes 
specific choices for all three of these, which give a par- 
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FIG. 1: Stationary Schwarzschild solution, Nester gauge. Shown is the region —0.16 < r < 4.0. The outer boundary is located 
at r = 9.84. 




FIG. 2; Stationary Schwarzschild solution, Lorentz gauge. 
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ticularly simple symmetrizable hyperbolic system with 
all propagation either along the light cone or along the 
tetrad congruence. All we do in this paper is adapt 
the variables to explicit spherical symmetry, and com- 
pare two alternative tetrad gauge conditions, the Nester 
and Lorentz gauges. The fairly trivial choices of vari- 
ables that we have explored {Bt = 1/i? versus R and 
Bn = versus A) have some effect on accuracy, but 
no real effect on stability. 

We have discovered potential problems with both 
Nester and Lorentz gauges which cannot be fixed by play- 
ing with the equation formulation. The Lorentz gauge is 
not consistent with stationary solutions of the evolution 
equations which have asymptotically hyperbolic spacelike 
hypersurfaces. With the Nester gauge the tetrad congru- 
ence can develop a kind of gauge shock, with steep gra- 
dients in the acceleration, for certain initial conditions. 
These issues would remain even if other problems could 
be fixed with a change in equation formulation. 

In experiments with various boundary conditions on 
the evolution equations, we found that it is critically im- 
portant to use some form of constraint-preserving bound- 
ary conditions. Otherwise constraint errors are gener- 
ated at the boundary and propagate into the grid at the 
speed of light. The particular scheme we implement is 
not as elegant as some that have been proposed in the 
literature, but seems to work well, at least in this sim- 
ple one-dimensional context. The other key to bound- 
ary conditions is to arrange to have the " zero- velocity" 
modes along the tetrad congruence propagate out of the 
grid at the boundary. This is no problem at the inner 
excision boundary, but requires expansion of the tetrad 
congruence at the outer boundary. Since we keep our 
constant-time hypersurfaces close to being orthogonal 
to the tetrad congruence, this is equivalent to having 
asymptotically hyperbolic hypersurfaces bending toward 
the future. With our choice of sign for extrinsic cur- 
vature, this means positive extrinsic curvature for the 
hypersurface. 

For the Nester gauge, there is only a rather nar- 
row range of initial conditions which give rise to stable 
and reasonably accurate solutions of the evolution equa- 
tions. For stationary initial conditions, if the parameter 
Uq = R Kt is too close to zero, the acceleration and the 
radial extrinsic curvature have very large gradients 
near the horizon, which is bad for stability. However, Uq 
needs to be much less negative than —0.25 in order that 



the hyperbolic character of the hypersurface be strong 
enough to give stability. There is a rather narrow window 
around Uq — —0.14 where constraint and other errors do 
not grow substantially. We have not found any initial 
conditions for which the evolution in the Lorentz gauge 
is stable. 

The importance of positive extrinsic curvature for sta- 
bility has also been noted in a much more rigorous ana- 
lytic treatment of the stability of the Weyl tensor evolu- 
tion equations by Frauendiener and Vogel [Zlj . 

Finally, we note that in the context of black hole evolu- 
tion it is critical to control the relation of the inner edge 
of the coordinate grid to the apparent event horizon, and 
highly desirable to keep the outer edge of the grid at a 
constant physical radius. We do this through the bound- 
ary conditions on our minimal deformation equation for 
the shift vector. Failure to do this leads to severe numer- 
ical problems. 

We conclude that the WEBB equation formulation 
as presented in jl5|, with either the Nester or Lorentz 
gauge conditions, does not seem particularly promising 
for black hole evolutions. While the Nester gauge can 
give good results for some initial conditions, these re- 
sults are not very robust. The Lorentz gauge does not 
seem to be viable at all. However, it may be possible, by 
adding constraint damping terms to the evolution equa- 
tions, to improve accuracy and stability for a wider range 
of initial conditions, while preserving a reasonable hyper- 
bolic structure. We have begun to explore other equation 
formulations and gauge conditions in the context of the 
tetrad framework, some of which seem considerably more 
promising. 
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FIG. 5: Convergence of energy constraint errors ai t — 40 M for stationary IVP and Nester gauge. The results at resolution 
dr = 0.01 have been muhiplied by 4, and those at dr = 0.005 have been muhiplied by 16. dr/dt = 10.0 for all resolutions. The 
lapse and shift are reset every dt — 0.002. 




FIG. 6: CMC initial conditions and Nester gauge: With Ko — 0.4, evolves pathologically, but with Ko = 0.2, approaches 
a stationary solution (shown is region near event horizon). *These times are for r = 9.84. 



